BCVEGPY2.0: A upgrade version of the generator BCVEGPY 
with an addendum about hadroproduction of the P-wave B c states 

Chao-Hsi Chang 1 ' 2 *, Jian-Xiong Wang 3 t and Xing-Gang Wu 3 * 
1 CCAST (World Laboratory), P.O.Box 8730, Beijing 100080, P.R. China} 
2 Institute of Theoretical Physics, Chinese Academy of Sciences, 

P.O.Box 2735, Beijing 100080, P.R. China. 
^Institute of High Energy Physics, Chinese Academy of Sciences, 
P.O.Box 918(4), Beijing 100049, P.R. China. 

Abstract 

The generator BCVEGPY is upgraded by adding the hadroproduction of the P-wave excited B c 
states (denoted by B* JL=1 or by hs c and xb c ) and by improving some features of the original one as 
well. We denote it as BCVEGPY2.0. The P-wave production is also calculated by taking only the 
dominant gluon-gluon fusion mechanism (with the subprocess gg — * B* JL=1 +c+b being dominated) 
into account as that for 5-wave. In order to make the addendum piece of the upgraded generator 
as compact as possible so as to increase its efficiency, we manipulate the amplitude as compact 
as possible with FDC (a software for generating Feynman diagrams and the algebra amplitudes, 
and for manipulating algebra formulae analytically etc) and certain simplification techniques. The 
correctness of the program is tested by checking the gauge invariance of the amplitude and by 
comparing the numerical results with the existent ones in the literature carefully. 
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PROGRAM SUMMARY 



Title of program : BCVEGPY 
Version: 2.0 (December, 2004) 

Program obtained from : CPC Program Library or Institute of Theoretical Physics, Chi- 
nese Academy of Sciences, Beijing, P.R. China: www .itp.ac.cnj~ zhangzx / BCV EG PY2.0 / . 

Reference to original program : BCVEGPY1.0 

Reference in CPC : Comput. Phys. Commun. 159, 192(2004) 

Does the new version supersede the old program : Yes 

Computer. Any computer with FORTRAN 77 compiler. The program has been tested 
on HP-SC45 Sigma-X parallel computer, Linux PCs and Windows PCs with Virsual Fortran. 

Operating systems : UNIX, Linux and Windows. 

Programming language used : FORTRAN 77/90. 

Memory required to execute with typical data : About 2.0 MB. 

No. of bytes in distributed program, (including PYTHIA6.2) : About 1.1 MB. 

Distribution format : Compressed tar file. 

Keywords : Event generator for hadronic production, B c meson, P-wave B c states, 
Upgraded version 

Nature of physical problem : Hadronic production of B c meson and its P-wave excited states. 
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Method of solution : The code with option can generate weighted and un-weighted events. 
For jet hadronization, an interface to PYTHIA is provided. 

Restrictions on the complexity of the problem : The hadronic production of cfr-quarkonium 
in S"-wave and P-wave states via the mechanism of gluon-gluon fusion are given by the 
so-called 'complete calculation' approach. The less important contributions from the other 
mechanisms have not been included. 



Typical running time : It depends on which option one runs to match PYTHIA when gen- 
erating the B c events. Typically, For the hadronic production of the S-w&ve cfr-quarkonium, 
if IDWTUP=1, then it takes about 20 hour on a 1.8 GHz Intel P4-processor machine to 
generate 1000 events; however if IDWTUP=3, to generate 10 6 events, it takes about 40 
minutes only. For the hadronic production of the P-wave cfe-quarkonium, the time will be 
almost two times longer than the case of the S'-wave quarkonium. 



LONG WRITE-UP 



I. INTRODUCTION 

Recently, we have completed the hadronic production of the P-wave states of B c 
that has not been included in the existent version BCVEGPY1.0 3] yet, so we upgrade the 
generator here. The studies in Refs. involve not only the contributions from the color- 
singlet components |(c6)i( 1 Pi)) and |(c6)i( 3 Pj)) (J = 1,2,3) but also those from the color- 
octet components | (c&)s( 1 »S'o)(7) and \(cb)s( 3 Si)g) [4] to the P-wave production. Especially, 
it is found that the P-wave production at Tevatron and LHC can be so big even as that in 
the same order of the magnitude for the S'-wave production in certain kinematics regions. 
Therefore, to compare with the original generator, in addition to some improvements on the 
program for the 5-wave B c (and B*) production, the main fresh feature of the upgraded 
version BCVEGPY2.0 is that it can generate the hadropro duct ion of the P-wave excited 
B c states (denoted by B* JL=1 or by /ie c , x% c i Xb c an d XbJ as wen - At such high energies 
of Tevatron and LHC, the gluon-gluon fusion mechanism i.e. that via the subprocess gg — > 
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B*j L=l + c + b, is also dominant for the P-wave production , hence in the addendum 
we only take the gluon-gluon fusion mechanism into account. 

Therefore, in fact, here we just need to rewrite the program used in Refs. into the 
format according to the environment PYTHIA 5] as the original version BCVEGPY1.0 
properly. 

Besides that the hadronic production of the P-wave B c states is added, for convenience 
and increasing efficiency, some improvements in BCVEGPY2.0 are made. They can be 
summarized as: i). an improved way to record the B c events is available, i.e., all the 
possible information, such as various distributions of the B c events with proper kinematic 
cuts, can be recorded just by running the program once; ii). an improved way is available 
to generate the useful sampling importance function just by running VEGAS program [f| 
carefully; iii). an 'extra switch' is added, i.e. whether one switches on or off to adopt the 
newly improved parton distribution functions (PDFs) (other than the inner PDFs provided 



y the latest version of the PDFs 
|, MRS conveniently; iv). the 



by PYTHIA itself) to the simulation, so that one can apt 
provided by the several groups such as CTEQ 7], GRV |j 
color flow of the processes is re- written (as an important issue for the new version), so as to 
well-match that in PYTHIA. 

In BCVEGPY1.0, the color flow is treated according to the usual fundamental- 
nation deposition ,10]. Since the co.or flow for a given hard scattering snb- 
process in PYTHIA is treated in the way as that in Ref.^jJ, i.e., the color fundamental- 
representation decomposition is replaced by the color-flow decomposition equivalently under 
the large N c limit. Such a color-flow decomposition has a more intuitive physical interpre- 
tation and practically it is easier to deal with than the usual fundamental-representation 
decomposition. Hence in order to well-match with PYTHIA and to be easy to deal with, 
especially, with the shower Monte Carlo for QCD jets, the replacement to treat the color 
flow as that in Ref. [ill] for the relevant sub-processes is adopted in BCVEGPY2.0. 

In order to increase the efficiency and to make the P-wave production amplitude as 
compact as possible, we manipulate the amplitude with the help of some techniques and the 
software FDC (a software for generating Feynman diagrams and the algebra amplitudes, 
and for manipulating algebra formulae analytically etc jl2j). To guarantee the correctness 
of the program, as Refs.pl Q , we have made numerical checks on the gauge invariance of 
the amplitude and numerical comparisons between the present obtained results and those 
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in the literature carefully. 

The paper is organized: following Introduction, in Section II we outline the new features 
of BCVEGPY2.0 in detail, and explain how to use and test (check) the program. In Sec. Ill, 
we present some discussions on the generator and make a summary. 



II. THE PROGRAM: BCVEGPY2.0 



Since in hadron collisions at high energies, the gluon-gluon fusion mechanism for the 
hadronic production of the B c meson is dominant over the other mechanisms, thus in 
BCVEGPY1.0 [3], the hadronic production of the S- wave B c and B* mesons via the gluon- 
gluon fusion mechanism has been accomplished. While the present target for the new version 
BCVEGPY2.0 now essentially is to add the hadronic production of the P-wave B c states 
to the generator BCVEGPY. In BCVEGPY2.0, for the production of the four P-wave 
B* JL=1 (Iib c , X°b c i Xb c an d XbJ states, not only the contributions from color-singlet com- 
ponents |(c6)i( 1 Pi)) and |(c6)i( 3 Pj=o,i,2)) but also those from the color-octet components 
\(cb)s( 1 So)g) and |(c6) 8 ( 3 S , i)^) are included. The formulation and additional techniques in 
calculating the hadronic production of P-wave states have been given in Refs. ]} m with 
explanations, hence for simplicity, we do not repeat them here. The readers may refer to 
Refs.Q, 3| for details if they are interested in them. 

To supersede the old version BCVEGPY1.0 by the new version BCVEGPY2.0 conve- 
niently, BCVEGPY2.0 is also written into a Fortran package, and it can generate the S"-wave 
and the P-wave B c states accordingly with proper options. In BCVEGPY2.0, for compar- 
ison, the hadronic production of S'-wave states P c ( 1 S'o) and B*( 3 Si) via the less important 
light quark-antiquark annihilation subprocess, q + q —>■ B C (B*) + c + b, (the annihilation 
mechanism) has also been included. Since BCVEGPY2.0 is also written in the format as 
PYTHIA (including common block variables), thus all the functions of PYTHIA can be 
utilized in connection with the use of BCVEGPY2.0. 



A. New features in BCVEGPY2.0 



BCVEGPY2.0 is written in the same structure as that of BCVEGPY1.0, with some 
improvements in addition to the P-wave production being merged. It contains the following 
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new features: 

• The amplitudes for the hadronic production of the color-singlet corresponding to the 
four P-wave states, B* JL=1 or 1 Pi and 3 Pj (J=0,l,2), are included. Four files to calculate 
the amplitudes for the P-wave states are added: plplamp.for, pjOamp.for, pj lamp. for and 
pj2amp.for, which are for the amplitudes of 1 Pi and 3 Pj (J=0,l,2) respectively. While the 
necessary subroutines and functions for calculating the square of the amplitudes for all the 
P-wave states are put into the files psamp.for and lorentz.for. In fact, the file lorentz.for 
includes basic subroutines, which are to calculate the polarization vectors and the polariza- 
tion tensor of the heavy quarkonium, and to generate short notations for typical expressions 
of useful contracts of the vectors and/or tensors so as to make the program more compact. 
Three extra head files: inclamp.f, inclcon.f and inclppp.f are included in the package, which 
contain the frequently used common blocks for calculating each hard scattering amplitude. 
Note here, some new techniques to compute the P-wave amplitude j]J, such as to expand 
the amplitude by 'independent and elementary fermion strings' analytically with the help 



of the FDC program 



12j etc., are used instead of the original helicity ones that are used in 



computing the S'-wave amplitude. 

• To give full estimates of the P-wave production up to v 2 , we also take into account 
the contributions to the hadronic P-wave production from the color-octet components 
\(cb)s{ 1 S )g) and |(c6) 8 ( 3 Si)g). In the cases of the hadronic production gg — > (cfo)i + b + c 
with (c6)i in color-singlet, |(c6)i( 1 S'o)) or |(c6)i( 3 Sj)}, there are only three independent 
color factors (i.e. three independent color flows UJ$). While in the cases of the P- 
wave hadronic production due to the color-octet components, such as gg — > (c6) 8 + b + c 
with (cb) s m color-octet, there are totally ten independent color factors 2]. We shall show 
in more detail about how to deal with the color-flows in the next subsection. The hard 
scattering amplitudes corresponding to the color-octet components can be directly read 
from the Eqs.(4-8) in Ref.0] but the color-singlet matrix elements (or the wave function 
at origin) are replaced by the color-octet matrix elements i^)\x\T d ip c (a i H aH)'4 } lT d Xb\fy and 
{0\xWT d ip c {a^aH^la^XblO) accordingly. We can approximately relate the color-octet 
matrix elements to the color-singlet one by NRQCD scaling rules |4| as the method de- 
scribed in Ref . |l^ . Applying the heavy-quark spin symmetry, the above two production 
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matrix elements are related by Q] 

(0\xWT d i> c (aW)4^T d Xb\0) = (2 J + 1) • (0| X jT<ty c ^ H a H )^ c T d Xb |0)[1 + 0(t; 2 )], (1) 

where J is the total angular momentum of the hadron state. More specifically, based on the 
velocity scale rule [4j, we estimate 

(0\xlT d i> c (a? H a H )il>lT d Xb\0) ^ A s (v) 2 ■ {^ h ip c ^ H a H )^ c x b \^) 

~ A^^.^Olxfel^o))! 2 .^^ 4 )] (2) 

where the second equation comes from the vacuum-saturation approximation, v is the 
relative velocity between the constitute quarks in the bound states. A$(v) is of the order 
v 2 or so, and one can take it to be within the region of 0.10 to 0.30, which is in consistent 
with the identification: A$(v) ~ a s (Mtj). 

• The 5-wave ( 1 Sq and 3 Si) hadronic production via the less important mechanism, i.e. 
the light quark- ant iquark annihilation mechanism with the subprocess q+q — > B c (B*)+c+b, 
is also included in BCVEGPY2.0. In the light quark-antiquark annihilation mechanism, the 
technique as that in the gluon-gluon fusion one [jj is applied, but it is much more simple 
(only seven Feynman diagrams are needed to be considered), so we do not present any detail 
for the calculation here. Numerical results of the calculation show that its contributions to 
the production are smaller than those from the gluon-gluon fusion (~ 1%) 



15. 



16]. 



• For convenience, in order to record more dcitci £it one run which may interest us, we add 
24 extra files to record the information of the generated events, and then all these output 
files are put in a subdirectory named data. For example, the generated files, pt005y.dat, 
pt020y.dat, pt035y.dat, pt050y.dat and ptl00y.dat are used to record the rapidity distribu- 
tions with different p t cuts, i.e. p tcu t = 5, 20, 35, 50, lOOGeV, respectively. The user may 
record other interested information in a similar way. 

• With the option IVEGASOPEN=l, one can use the VEGAS program j(| to achieve an 
important sampling-function, which may be used to increase the efficiency of the Monte Carlo 
simulation j^]. Otherwise one may set IVEGASOPEN=0. Furthermore, in order to improve 
the precision of the existed important sampling-function obtained by previous VEGAS run, 
we add one more input parameter IVEGGRADE, i.e., with IVEGGRADE=1, one can use 
the existed sampling important function recorded in the file grade.dat to generate a more 
precise important sampling-function and record it in a new file newgrade.dat. Otherwise 
one should set IVEGGRADE=0. 
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• PYTHIA itself provides sixteen parton distribution functions (PDFs). It is known that, 
in fact, the PDFs are obtained through global fitting of the experimental data and then in- 
voluted to the requested characteristic scale Q 2 by a standard way of perturbative Quantum 
Chromodynamics (pQCD). Several groups, CTEQ 0, GRV 8( and MRS 9] etc, devote 
themselves to offer accurate PDFs to the world and upgrade them continually with fresh ex- 
perimental data being available. BCVEGPY2.0 provides a switch (IOUTPDF) to determine 
which type of PDF one uses: if (IOUTPDF=l), then one can use one of the latest version of 
PDFs from the above three groups with a proper value for the parameter (IPDFNUM) to do 
the hadronic production; otherwise, with (IOUTPDF=0), one can directly use a certain exis- 
tent PDF (by setting the value of the PYTHIA parameter MSTP(51)) provided in PYTHIA 
to do the hadronic production. For convenience, we have downloaded all the source files 
for the latest version of the leading order (LO) PDFs of those three groups, i.e. CTEQ6L, 
GRV98L and MRST2001L, and saved them into a single fortran file outerpdf.for. The three 
corresponding data files (cteq61.tbl, grv981o.grid and lo2002.dat), which are necessary when 
running the file outerpdf.for for these three PDFs, have been also added in the package. In 
fact, it is shown by numerical calculations that the uncertainties for the estimates due to the 
different PDFs can reach up to the level ~ 10% in both cases for the S'-wave and P-wave 

nn 

hadronic production |ll 115]. 



• For convenience, we have written an additional file parameter. for to set the initial 
values of the parameters shown in TAB. I. One may fix the value of 'IREADDATA' in 
parameter. for as zero ('IREADDATA=0') and then one returns back the original way 
(that in version 1.0) to set the input parameters by reading the parameter values from the 
date file: totput.dat. In practice, if one want to generate a large number of results with 
various input parameters, we would like to suggest him/her to adopt the way of reading the 
input data directly from the date files accordingly. In this way, only the input files need to 
be changed and all the source files (*.for) need not to be changed, so we can save a lot of 
time on making the source files and furthermore, all the works can be written in a single 
executable file and be executed automatically. 

In summary, in the new package BCVEGPY2.0, there are twenty one files totally, 
which includes thirteen main source files: bcvegpy.for, foursets.for, genevnt.for, lorentz.for, 
outerpdf.for, plplamp.for, parameter. for, pjOamp.for, pj lamp. for, pj2amp.for, psamp.for, 
py6208.for (PYTHIA version 6.208), ssamp.for; three date files for the corresponding three 
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TABLE I: The parameter values in the sequential order in the totput.dat file. 
PMBC PMB PMC FBC 

PTCUT ETACUT ECM IBCSTATE IGENERATE IVEGASOPEN 
NUMBER NITMX 
NQ2 NPDFU NEV 

ISHOWER MSTP(51) IDWTUP MSTU(lll) PARU(lll) 
ISUBONLY SUBENERGY IGRADE 
INUMEVNT IVEGGRADE IQQBAR IQCODE 
IOUTPDF IPDFNUM IOCTET COEOCT 



PDFs: lo2002.dat, cteq61.tbl, grv981o.grid; three head files: inclamp.f, inclcon.f, inclppp.f; 
one input data file: totput.dat and finally a simple help file: readme.dat. 

With all these new features in mind, the sequential order and the format of the parameters 
in the input data file, totput.dat, now changes to Tabled. Some new parameters specified in 
the input file are: 

• FBC=: is different from the meaning in BCVEGPY1.0, and now it's value is the radial 
wave function R(0) for S'-wave production or the first derivative of the radial wave function 
R'(0) for P-wave production. In our calculations, we take FBC=1.24GeV 3 / 2 , which cor- 
responds to ^(0) = 1/V^R(0) = 0.35Ge\/ 3 / 2 for S-wave, and FBC=0.448Ge\/ 5 / 2 , which 
corresponds to ^'(0) = ^3/4^(0) = 0.219GeV 5 / 2 for P-wave. 

• IBCSTATE=: indicates the S- and P-states of B c meson: IBCSTATE=1, for P c [%]; 
IBCSTATE=2, for P c *[ 3 Si]; IBCSTATE=3, 4, • • • , 6, for the four P-wave states P^Pi], 
P*[ 3 P ], P*[ 3 Pi] and B*[ 3 P 2 ] respectively. 

• INUMEVNT=: indicates whether to switch on or off to keep the event number distri- 
butions other than the differential distributions of p T , y and etc.. INUMEVNT=1 means to 
switch on and INUMEVNT=0 means to switch off. 

• IVEGGRADE=: indicates whether to switch on or off to use the existed grade (impor- 
tance sampling function) before running VEGAS. IVEGGRADE=1 means to switch on and 
IVEGGRADE=0 means to switch off. With this parameter, one can improve the efficiency 
of the importance sampling function. 

• IQQBAR=: indicates whether to switch on or off to use the less important light quark- 
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antiquark annihilation subprocess q + q — > (c6) + c + 6 to generate the (cfe)-quarkonium 
events, where (c6)-quarkonium is in 1 S or 3 Si state. IQQBAR=1 means to switch on and 
IQQBAR=0 means to switch off. When 1QQBAR=0, it will use the default gluon-gluon 
fusion subprocess to generate the (c6)-quarkonium events. 

• IQCODE=: indicates the species of the light quark for the light quark-antiquark annihi- 
lation mechanism. Under the option IQQBAR=1, IQCODE=l means u and u, IQCODE=2 
means d and d, while IQCODE=3 means s and s used in the quark-antiquark annihilation 
subprocess. 

• IOUTPDF=: indicates whether to switch on or off to use the three latest version of PDFs: 
CTEQ6L, GRV98L and MRST2001L offered in the program for the hadronic production. 
IOUTPDF=l means to switch on, while IOUTPDF=0 means to switch off and then the 
inner PDFs of PYTHIA is used to generate the events. 

• IPDFNUM=: indicates which one of the three latest version of PDFs: CTEQ6L, 
GRV98L and MRST2001L is used for the hadronic production. It comes into operation only 
under the option IOUTPDF=l. IPDFNUM=100 means to use GRV98L; IPDFNUM=200 
means to use MRST2001L; IPDFNUM=300 means to use CTEQ6L. 

• IOCTET=: indicates whether to switch on or off to use the subprocess gg — > (cb)$ + b + c 
to generate the color-octet components \(cb)s( 1 So)g) and | (c6)s( 3 *S'i) (?) . IOCTET=l means 
to switch on and IOCTET=0 means to switch off. 

• COEOCT=: the value for A s {v), where A s {v) is defined in Eq.©. 



B. On the color flow for the hardronic production of B c meson 

A string-based fragmentation scheme for the cross sections, such as that in the Lund 
model and used in PYTHIA, is needed for different color flows. In order to merge well 
with the shower Monte Carlo used in PYTHIA, in BCVEGPY2.0 we adopt the color-flow 



11] to analyze the production. In Ref. 



Ill, 



decomposition that has been introduced in Ref. 
it is demonstrated that the color-ordered amplitudes (i.e. amplitudes pertain to the cor- 
responding color-flows), that appear in the color-flow decomposition, are identical to those 
appear in the usual fundamental-representation decomposition under the large N c limit. So 
practically, one can first use the color-flow decomposition to analyze how many independent 
color-flows we need for a concerned subprocess and then transform the results into the usual 
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FIG. 1: Color flow diagrams for the color-singlet production based on the color-flow 



decomposition 111]. Each pair of indices ik and jj~ corresponds to an external gluon, i.e. k = 1 
is for gluon-1 and k = 2 for gluon-2. i and j are the decomposed color indices for the outgoing c 
and b respectively. 



fundamental-representation decomposition according to the color-flow Feynman rules. In 
this way, we can read out the color-ordered amplitudes easily from our original results in 
Refs.QjJEJ. 

Since the less important light quark-antiquark annihilation mechanics with the subprocess 
q + q — > B C (B*) +c + b has only two independent orthogonal color-flows, so we do not repeat 
the results here. We show the color-flows for the dominant gluon-gluon fusion subprocess in 
detail as the follows. In BCVEGPY2.0, we need to deal with the color-flows for the dominant 
gluon-gluon fusion subprocess gg — > (cb) + b + c in two cases: one is for (c5)-quarkonium 
in a color-singlet state (i.e. B C ,B* and B* JL=1 ) and the other is for (c6)-quarkonium in a 
color-octet state (i.e. \(cb) 8 ( 1 S )g) and \(cb) s ( 3 Si)g)). 

• The first case: (c6)-quarkonium in a color-singlet state. In this case, we does not need to 
consider the color structure of the (cfe)-quarkonium and one may find that there are only 
three types of independent color flows as shown in Fig^J 



(3) 



where (i n ,jn) (n=l,2) are the color indices for the gluon-1 and gluon-2 respectively, i and 
j are color indices of the outgoing c and b. In the conventions, the lower indices i and i n 
(n=l,2) transform under the anti fundamental-representation of SU(N C ) and upper indices 
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FIG. 2: Color flow diagrams for the color-octet production based on the color-flow 



decomposition 111]. Each pair of indices ik and jk corresponds to an external gluon, i.e. k = 1 
is for gluon-1, k = 2 for gluon-2 and k = 3 for the color-octet (c6)-quarkonium. i and j are the de- 
composed color indices for the outgoing c and b respectively. The diagrams for c n+ \ [n = 1, 2, • • • , 5) 
can be directly obtained by gluon exchange. 

j and j n (n=l,2) under the fundamental-representation. Note here that the color bases are 
not normalized as done in Ref . , because we only need the relative probabilities among 
each color flow. According to the color-flow Feynman rules [llj], the color-flow decomposition 
can be related to the fundamental-representation through the following transformation: 



(T a T%, c 2 



(8 ij Tr[T a T b }). 



(4) 



Equivalently, the total amplitude of the process can be written in the fundamental- 
representation decomposition in the following way: 



M = {T a T b ) ij M 1 + {T b T a ) ij M 2 + (5 ij Tr[T a T b })M 3 , 



(5) 



where T a and T b are the color matrices of the gluon-1 and gluon-2; i,j = 1,2,3 are color 
indices of the outgoing quarks c and b respectively. Mi are the color-ordered amplitudes 



HQ. 



respectively, which can be directly read out from Refs. 
• The second case: (cfe)-quarkonium in the color-octet state. In this case, we need to consider 
the color structure of the (c6)-quarkonium itself, and one may find that there are ten types 
of independent color flows totally as shown in Fi^21 
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cs = (6Ml6i>5i), c 9 = {St <W4) > no (6) 

where 13 and j'3 are the decomposed color indices of c and b in the (c&)-quarkonium. Ac- 
cording to the color-flow Feynman rules the color-flow decomposition can be related to 
the fundamental-representation through the following 'correspondences': 

Ci ^ (T b T a T d )ij, c 2 -> (T a T b T d )ij, c 3 -> (^Trp-T^), 
c 4 -> (^Tr[T 6 T d ]), c 5 -> (^Tr[T 6 T a T rf ]), c 6 -> (^Tr[T a T 6 T rf ]), c 7 -> {T d T b T% v 

c 8 ^(T d T a T%, c 9 ^(T a T d T%, c 10 -> (T b T d T a ) l]} (7) 

where T d stands for the color matrix of the color-octet (c6)-quarkonium. Equivalently, the 
total amplitude of the process can be written in the fundamental-representation decompo- 
sition in the following, 

M = (T h T a T d )ijMi + (T a T b T%M 2 + (T b Tr[T a T d })M 3 + (T l a j Tr[T b T d })M 4 
+ (5ijTr [T b T a T d ])M 5 + (5 i:j Tr[T a T b T d })M6 + (T d T b T a )ijM 7 + (T d T a T b ) i:i M 8 
+ (T a T d T b )ijM 9 + (T b T d T a )ijM 10 . (8) 

Mi are the corresponding color-ordered amplitudes, which can be directly read out from the 
results in Refs. 

• In both cases, when squaring the amplitude and summing over colors, the leading terms 
in 1/N C obtained by the square of each color flow. The cross terms between different color 
flows are suppressed by a certain power of N c comparing with the leading terms and under 
the large N c limit they can be safely neglected. Therefore, we can obtain the probability 
for each color flow, which are just |Mj| 2 accordingly, with i=(l,2,3) in color-singlet case and 
i=(l,2,- • -,10) in color-octet case. 



C. Checks for the generator 

The amplitude for S'-wave production has been checked well in the original version already, 
thus the check for the whole Fortran package BCVEGPY2.0 essentially is on the new contents 
in comparison to the old version. As for BCVEGPY2.0, it is checked by examining the 
gauge invariance of the amplitude for each production, S'-wave and P-wave of the B c meson, 
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TABLE II: The total cross section for the hard subprocess g + g — > B* JL=1 + b + c (gluon fusion 
into P-wave excited states B* JL=1 ) at different CM. energies. The input parameters are taken as 

n 

those used in Ref.|17j|: m& = b.OGeV, m c = 1.7GeV, M = 6.7GeV and the running a s is fixed to 
0.2 etc. 



CM. energy (GeV) 


20GeV 


40GeV 


60GeV 


80GeV 


lOOGeV 


200GeV 


a( 1 P 1 )(pb) 


0.367 


0.743 


0.657 


0.538 


0.439 


0.195 


a( 3 P )(pb) 


0.184 


0.207 


0.175 


0.141 


0.114 


0.0496 


a( i P 1 )(pb) 


0.346 


0.598 


0.503 


0.402 


0.324 


0.139 


a( 3 P 2 )(pb) 


0.721 


1.49 


1.31 


1.06 


0.862 


0.374 



10° 




10 I 1 1 1 1 1 1 1 1 1 ' i i 

5 10 15 20 25 30 35 40 45 50 1 1 1 1 1 1 ! 1 

//-\ v t\ -3-2-10123 

P t (GeV) y 

FIG. 3: pt and y differential distributions for the subprocess, where cb is in P-wave states. The 
dotted line, the dash-dot line, the dash line and the solid line are for 1 P\, 3 -Po> 3 Pi and 3 P2, 
respectively. 

individually, i.e. the amplitude vanishes when the polarization vector of an initial gluon 
is substituted by the momentum vector of this gluon. Numerically we find that the gauge 
invariance is guaranteed at the computer ability (double precision). 

For further checks, we have performed several test runs. With the option ISUBONLY=l, 
we obtained the transverse momentum p t and rapidity y distributions and the total cross- 
sections for the sub-process gg — > B*j L=1 + c + b. The results by taking the same parameters 
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TABLE III: Generation parameters used in the sample generation. 



INITIAL PARAMETERS .. 

b 



GET THE RESULT FOR Be IN 3 R 



Be IN COLOR-SINGLET STATE 

GENERATE EVNTS 30000000 FOR TEVA ENERGY(GEV) 0.196£ + 4 

* USING SUBPROCESS: g + g^Bc + b + c * 

M_{Bc}=6.700 M_{B}= 5.000 M_{C}=1.700 f_{Bc}=0.2879 
Q2 TYPE= 3 ALPHAS ORDER= LO 

THE NEW VERSION OF OUTER PDFs WHICH ARE FROM WWW 
PDF: CTEQ6L; ALPHA IN LO 
USING PYTHIA IDWTUP= 3 
PTCUT =0.000 GeV 
NO RAPIDITY CUT 

USING VEGAS: NUMBER IN EACH ITERATION= 200000 ITERATION= 30 
END OF INITIALIZATION 



as those used in Refs.|l7l. Il8l| are shown in the TABLE|H]and in FIGH3 One may see that 
the results for the total cross-section of the hard subprocess g + g — »■ B* JL=1 + c + b with 
various CM. energies are agree with those in Ref. 171 (TABLE. I there) within the MC error, 

nn 

and the r>t and y distributions are consistent well with the ones obtained in Refs. 17, 18 . 

nn 

Since our present method is quite different from that in Refs. [13, I18L th e comparison is 
really a good check for the generator and also for the results in Refs. fl?!. In fact, in 
Refs. (17! Il {J], the derivative over the relative momentum q of the constitute quarks in the 
B c meson, which is necessary for P-wave production, is done numerically; while in our case, 
it is done analytically with the help of the FDC program |l2j |. 

Finally, we show a test run for the hadronic production of ( 3 -Po) cfe-quarkonium at the 
TEVATRON. When running the program, the initialization is shown as a screen snap-shot 
in TABLE ED Here the value of f Bc = 0.2879 corresponds to R'(0) = O.U8GeV 5 / 2 . The 
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0.04 



0.1 i 




20 40 60 80 100 -4 -2 2 4 



p t (GeV) y 



FIG. 4: Test run results of the hadronic production of 3 Pq state at TEVATRON. The transverse 
momentum pt and rapidity y distributions for the CTEQ6L by using a s in the leading order (LO) 
and adopting the characteristic energy scale squared Q 2 = pf + M 2 c ^ . 



BEGIN 

X 



BASIC INPUT: NEV, NUMBER, ITMX.VEGASOPEN 



p^RnPFTJ>-^^H 4^ ITIALIZE GRADE 




IT=0; IT=IT+1 



NUM=0; NUM=NUM+1 



CALL PHPOINT() 
TOTFUN(X.WT) 




FIG. 5: Flow chart for using VEGAS. The meaning of the symbols can be found in the program. 



results are shown in FIG0J 



D. About VEGAS 



This subsection is devoted to give a supplementary explanation on VEGAS. 
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By using VEGAS, one can obtain an important sampling-function (or simply called grade 
function) into the file grade.dat, and as it is well known, by using such function to do the 
event simulation, the efficiency can be greatly increased |5[. A modification to the original 
version in Ref. 6( is done so as to record the obtained grade function into the data file 
grade.dat. Such small trick is useful, i.e. when one wants to generate the events with the 
same input parameters, then one can directly apply the grade function obtained by the 
previous running and does not need to run VEGAS again. The flow chart for this part is 
presented in FIG|SJ where as mentioned above, one can use the existent grade function to 
do the initialization further and then to obtain a more precise grade function. 

In calling VEGAS, we need to input two parameters, i.e. NUMBER (maximum total 
number of the times to call the integrand in each iteration) and ITMX (maximum number 
of allowed iterations). These two parameters should be adjusted in a suitable way so as to 
obtain a grade function in good precision and also to save the CPU as much as possible. 
One subtle point is that, in VEGAS, the default number of the column bin is 50 (i.e. the 
region of (0,1) is divided into 50 pieces); by taking a bigger proper value for the column 
bins, one can improve the production efficiency but need more iterations (or numbers) to 
obtain a stable result. In practice, the values for NUMBER, ITMX and the column bins 
should be carefully chosen to obtain a best important sampling-function within the least of 
time. In our program, we keep the column bins to be 50 (according to our experience, it is 
enough); however the user who wants to do some very precise studies, a proper variation of 
column bins might be needed. 



III. CONCLUSIONS 



An upgrade version of B c meson generator BCVEGPY2.0 is given. In this new version, 
the hadronic production of P-wave states up to v 2 , that include both the color-singlet com- 
ponents B* JL=1 and the color-octet components \ (cb)s( 1 So)g) and | (c6) 8 ( 3 S'i)g), is available. 
Since almost all of these P-wave excited states decay to the ground state B c meson through 
the electromagnetic transition (if the phase space for strong interaction transitions is not 
allowed) finally, thus in fact the P-wave production provides additional contributions to 
the hadronic production of the B c meson. In the new version, the hadronic production of 
.S-wave cfe-quarkonium from the quark- ant iquark annihilation subprocess is also included for 
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comparison, although the contributions from the quark-antiquark annihilation mechanism to 
the B c hadronic production, in comparison with those from gluon-gluon fusion mechanism, 
can be ignored safely at LHC/TEVATRON energies in the allowed kinematic region. 

To achieve a compact program, we apply the FDC program to manipulate the amplitude 
of the dominant gluon-gluon fusion subprocess for the production corresponding to the four 
color-singlet P-wave states The correctness of the program is guaranteed by various 
checks. The generator is interfaced with PYTHIA, which takes care of producing the full 
event and filling the standard PYTHIA event common block. In view of the prospects for 
B c physics at Tevatron and at LHC, the generator offers a platform for further experimental 
studies. 
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APPENDIX A: ADDITIONAL ROUTINES AND FUNCTIONS FOR CALCU- 
LATING THE HADRONIC PRODUCTION OF B* JL=1 

In this Appendix, additional subroutines and functions for calculating the amplitudes for 
the gluon-gluon fusion subprocess, gg — > B* JL=l + b + c, are explained. 
SUBROUTINE ampl_lpl(cc) 

Purpose: with 1=1, 2, • • • , 36 to compute the amplitude of the 36 Feynman diagrams of the 
gluon-gluon fusion subprocess correspondingly, with the (cfe)-quarkonium in X P\ state. 
Real*8 cc=: cc = ±^J 9R 2 ^ - 72 ^ j3 , where P'(0) 2 is the square of the first derivative of the 
radial wave function at the origin, and M is the mass of the (c6)-quarkonium. Whether it 
takes the plus or the minus sign is our inner convention for the 36 Feynman diagrams and 
it depends on how to decompose the Feynman diagrams. 
SUBROUTINE ampl_3p0(cc) 

Purpose: with 1=1, 2, ■ • • , 36, to compute the amplitude of the 36 Feynman diagrams of the 
gluon-gluon fusion subprocess correspondingly, with the (cfe)-quarkonium in 3 P state. 
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Real*8 cc=: cc = ^V ^^^iaMS ' wnere R'(0) 2 is the square of the first derivative of the 
radial wave function at the origin and M is the mass of the (c6)-quarkonium. Whether it 
takes the plus or the minus sign is our inner convention for the 36 Feynman diagrams and 
it depends on how to decompose the Feynman diagrams. 
SUBROUTINE ampl_3pl(cc) 

Purpose: with 1=1, 2, • • • , 36, to compute the amplitude of the 36 Feynman diagrams of the 
gluon-gluon fusion subprocess correspondingly, with (c6)-quarkonium in 3 P 1 state. 
Real*8 cc=: cc = ±\J 9R ' 2 ^ 2 - 14 ^ M 3 , where R'(0) 2 is the square of the first derivative of the 
radial wave function at the origin and M is the mass of the (c6)-quarkonium. Whether it 
takes the plus or the minus sign is our inner convention for the 36 Feynman diagrams and 
it depends on how to decompose the Feynman diagrams. 
SUBROUTINE ampl_3p2(cc) 

Purpose: with 1=1, 2, • • • , 36, to compute the amplitude of the 36 Feynman diagrams of the 
gluon-gluon fusion subprocess correspondingly, with (c6)-quarkonium in 3 P 2 state. 
Real*8 cc=: cc = ±\J 9R 2 ^ - 72 j^ 3 , where R'(0) 2 is the square of the first derivative of the 
radial wave function at the origin and M is the mass of the (c6)-quarkonium. Whether it 
takes the plus or the minus sign is our inner convention for the 36 Feynman diagrams and 
it depends on how to decompose the Feynman diagrams. 
FUNCTION amps20_lpl() 

Purpose: to compute the square of the whole amplitude for a particular combination of two 
gluons' polarization states and the polarization vector of l P\ cfr-quarkonium state. 
FUNCTION amps2_lpl() 

Purpose: to compute the square of the whole amplitude with the help of the FUNCTION 
amps20_lpl(), where the summation over the gluons' polarization states and the polarization 
vector of 1 P\ cfe-quarkonium state is done. 
FUNCTION amps20_3p0() 

Purpose: to compute the square of the whole amplitude for the production of 3 Po cb- 
quarkonium and for a particular combination of the two gluons' polarization states. 
FUNCTION amps2_3 P 0() 

Purpose: to compute the square of the whole amplitude with the help of the FUNCTION 
amps20_3p0(), where the summation over the gluons' polarization states is done. 
FUNCTION amps20_3pl() 
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Purpose: to compute the square of the whole amplitude for a particular combination of two 
gluons' polarization states and the polarization vector of 3 Pi cfr-quarkonium state. 
FUNCTION amps2_3pl() 

Purpose: to compute the square of the whole amplitude with the help of the FUNCTION 
amps20_3pl(), where the summation over the gluons' polarization states and the polarization 
vector of 3 Pi cfr-quarkonium state is done. 
FUNCTION amps20_3p2() 

Purpose: to compute the square of the whole amplitude for a particular combination of two 
gluons' polarization states and the polarization tensor of 3 P 2 cfe-quarkonium state. 
FUNCTION amps2_3p2() 

Purpose: to compute the square of the whole amplitude with the help of the FUNCTION 
amps20_3p2(), where the summation over the gluons' polarization states and the polarization 
vector of 3 P 2 c6-quarkonium state is done. 
SUBROUTINE genpolar3(p, ep, pm) 

Purpose: to generate the polarization vectors (ep) of a vector meson with mass pm and 
momentum p. 

real*8 p(4), ep(4,2:4): p(4) is the momentum for the vector meson, p(l) = E, p(2) = p x , 
p(3) = p y and p(4) = p z . ep(4,2:4) records three polarization vectors of the vector meson, 
each of which has four components. The explicit formulas of ep can be found in appendix 
A of Ref.Q. 

SUBROUTINE genpolar(p, ep) 

Purpose: to generate the polarization vectors (ep) of a gluon with momentum p. 
real*8 p(4), ep(4,2:3): ep(4,2:3) records two transverse polarization vectors of the gluon. 
SUBROUTINE gentensor(p, ep, pm) 

Purpose: to generate the polarization tensor (ep) for the meson in 3 P 2 state with mass pm 
and momentum p. 

real*8 p(4), ep(4,2:6): ep(4,2:6) records five polarization tensors of the meson in 3 P 2 state. 
The explicit formulas of ep can be found in appendix A of Ref . Q . 

There are also some accessory subroutines and functions that are used to simplify the 
amplitude and to make the amplitude more compact and not difficult to be understood, 
such as the SUBROUTINE genppp which is used to generate short notations for all the 
typical linear combination of dotted products of two momenta, thus for shortening the paper 
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we do not explain them precisely here. 
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